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Relaxation algorithm in description of 
superconducting structures. 
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Abstract 

Relaxation method is described on simple superconducting cases in one and two dimensions in Ginzburg-Landau 
(Gl) formalism. The structure of the algorithm is given for the case time independent and time dependent GL 
equation. The advantages and disadvantages of the algorithm are specified. Particular cases solved by the 
relaxation algorithm is unconventional Josephson junction (uJJ), modified uJJ, SQUID built on the base of uJJ 
and superconduting current limiter. The artifacts of relaxation algorithm are described. 
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Introduction to the relaxation algorithm 

Variational calculus is used in many branches of fundamental 
and applied science. In theoretical mechanics the minimiza¬ 
tion of action principle is used and it allows to derive the 
equation of motion. Maxwell equations can also be derived 
using variational derivative on action with respect to physi¬ 
cal quantities as electric or vector potential. Also quantum 
mechanics and theory of relativity can be formulated by mini¬ 
mization of action function and hence variational derivative 
with respect to desirable quantities will result in proper equa¬ 
tions of motion. Refined and theoretical review of relaxation 
method is given by Adler [?]. In this work we will consider 
only on effective formulation of superconductivity problems, 
which are similar in many ways to theory of superfluidity. We 
will concentrate on practical aspects of application of relax¬ 
ation method. In Ginzburg-Landau (GL) theory we have free 


energy functional and equations of motion are given by 


$F\Xjjjc J 
SXtjt,... 


( 1 ) 


where Xjjp is the physical field e.g. A-vector potential. 
More details on derivation of GL equations can be found in 
[?]. In numerical computations we adapt the scheme 


SF[XiJ,kJ _ AX,, U ,... 

SXijt,... rliJ ’ k ’- A,, ;X ... ’ 


( 2 ) 


where Aand is kept constant, while A 

is changing during the simulation. Using GL free energy 
functional F[y/,A] we obtain the Ginzburg-Landau equations 

8F[y,A\ _ ~ A yr(x) 8 F[xi/,A] ~ AA 

8y(x) m Ati ’ 8A 11 At 2 ' 


At first we consider the simple one dimensional form of 
Ginzburg-Landau equation of the form 

-—A x {x,t)) 2 -y^-)Y(x,t)= 0, (4) 

c at 

where y is some constant. Maxwell equation gives 


V x (V x A)) = Hoj 

curr (r) + jtt 0 £o 


dE(t) 

dt 


( 5 ) 


where j is the electric current density from GL theory. The 
superconductor-vacuum and superconductor-nonsuperconducor 
interface in yz plane is accounted in GL theory by condition 
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(is-= 0 or (f £ - = M*)* 

where b is some material depending on superconductor and 
non-superconducting material. For certain type of problems 
with translational symmetry in z direction and only A z com¬ 
ponent we can adapt the gauge VA(jc,y) = 0. Let us limit 
to one dimensional problems. Since we deal with numerical 
algorithm we obtain the approximate values of y/(x,t) across 
given lattice. Thus it will be helpful to introduce space error 
function for GL equation given as e(x) 

( v _ l(-s(fs-T A ^)) 2 + «W+^WIyWI 2 )yW)l 

C ' A; IV'MI + IV'mml 

( 6 ) 

where 11 (f m i n | is small constant. Since relaxtion method rely on 
many iterations with time it is also useful to introduce average 

y Nx e f x .\ 

error e av dehned as e av = I=1 N . It is quite straightforward 
to generalize e(x) and e av for two and three dimensions. 

Implementation of the relaxation 
algorithm 

Let us apply the relaxation method to the one dimensional 
GL equation 4, which is schematically depicted in Fig. 9. At 
first we make initial guess of vector potential A x (x) and 1 jf{x). 
We chose certain lattice of N x points and values of rji, rj 2 . 
At i and A i^. We compute initial electric current density j CU rr 
from GL theory and gradients of i \f and A x {x). Next step is 
the execution of N t times the loop: 

1. We compute j curr on every element of space lattice. 

2. We compute £yr(x,t), £zW(x,t), £A x (x,t) and ^A x (x,t). 
for each element of lattice. 

3. We compute the change of vector potential and y/ for every 
element of lattice 

A A x (x) = -Atiriij curr (*^) 5 

A yr(x) = At rj (a + J81 ip(x) | :2 + ^ (j lx ~ T A * W ) 2 ) W( x ) 



Figure 1. Distribution of superconducting order parameter 
(SCOP) y/(x) before and after simulation for the case of 1 
dimensional superconductor. 



Figure 2. Distribution of a coefficient. 


Various solutions obtained by the 
relaxation algorithm 

At first one dimensional GL equation given in 4 was solved 
with use of relaxation algorithm for the case of zero vector 
potential and a depicted in Fig. 2. The initial and final values 
of i f/ are given in Fig.l. The error relative error e(x) in the 
last step of simulation is depicted in Fig. 3. The final solution 
is obtained after the free energy F and average error e av are 
minimizing and get the saturation what is depicted in Fig. 4 
and in Fig. 5. Despite the fact that initial values of y/ 
were far from physical intuition proper numerical solution 
(corresponding to physical intuition) was obtained. 

Next example was solving two dimensional GL equation 
of the form 


4. We apply changes of Ay/ and AA x (x ) for each element of 
lattice 


f yf{x) -> y/(x)+Ayf(x), 

\ A z (x, y)-> A x (x) + A A x (x). 

5. We check the correctness of boundary conditions and make 
the adjustments in y/ so they are fulfilled. 

6. We make the correctness of certain physical constrains and 
make the necessary adjustments in y/ and A x held. 

7. We compute free energy F, numerical error e(x) and e av . 

The criteria necessary for obtaining the solution is the mini¬ 
mization and saturation of numerical error e(x), e av and free 
energy F. 


h 2 d 2 d 2 h 2 o Ae 2 

^ + ^ MX ' y) + { 2n, A ‘ (X '^ + 

a(x,y) + /3 (x,y) \ y/{x,y) | 2 ) y/{x,y) = 0 . 


( 7 ) 


At first we consider the case of zero magnetic held what 
implies A z (x,y) = 0. We continue the study on unconven¬ 
tional Josephson junction started by [?] when we place the 
non-superconducting element on the top of superconductor 
slab. This time the non-superconductor strip is placed in¬ 
side the superconductor. We also can consider placement of 
less superconducting bar on the more superconducting bar 
(more negative a coefficient). Such structures are depicted 
on the left and right side of Fig. 6, which are given by certain 
distribution of (x(x,y) held. 
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Figure 6. Distribution of oc(x,y) defining modified unconventional Josephson junction (left side) and half modified 
unconventional Josephson junction (right side). 




Figure 3. Distribution of relative error e(x) obtained in final 
step of simulation. 



Figure 4. Evolution of average error e av (t) with simulation 
steps. 


The presence of normal strip reduces superconducting 
order parameter (SCOP) in more significant way as in the case 
of weakening of the superconductor what is given on the left 
and right side of Fig. 7. The diffusion of Cooper pair from the 
superconducting strip into superconducting strip takes place. 
In the case with two superconductors superconducting order 
parameter diffuses from more superconducting region into 
less superconducting region. 

The situation with modeling superconducting mesoscopic 
structures becomes complicated when there is occurrence 
of non-zero vector potential and hence magnetic field. We 
can measure the error of y/ function as it was introduced by 
means of e(x) function. In order to measure the error of vector 
potential it is necessary to consider the equation for vector 
potential. We have 

V 2 A z (x,y) = -CA z (x,y)\\f/(x,y)\ 2 , (8) 

where C is constant quantity depending on fundamental phys¬ 
ical constants. Quite obviously from the last equation 


C(x,y) 


M x ,y) \v(x,y )\ 2 


(9) 


Free energy F 



Figure 5. Free energy F vs iteration t. 


one should obtain the constant value C. In real numerical 
simulation C values would change and be position dependent 
as C(i,j) for the case of discrete lattice. Therefore the best 
criteria is the minimization of quantity C\ is defined (seems 
to be defined) as 


y(\£C(k,l)\ + \^C(k,l)\) 
\ \C(k,l)\ + \C min \ 


where summation is conducted over the all points of lattice 
(. N X: N y ) and C m i n = 10 -3 for example. One could define 
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Figure 8. The introduction of stabilizing regions, which has the imposed \lf(x,y) = 0 (left side). In order to confirm the stability 
of solutions the numerical noise is introduced to the relaxation algorithm. It helps to achieve the global minima for F functional 
as depicted on the right picture. Such situation is analogical to the case of random dropping balls in random hills potential. 


many similar functions as C\ that should be minimized in 
the iterations of the relaxation algorithm. Quite much similar 
reasoning could be presented in case of 3 dimensions. All 
considerations are valid for the case of time independent and 
time dependent GL equations. 

Artifacts and limitations of the relaxation 
method 

Relaxation method is quite stable especially when iteration 
step At is small. For obvious reasons it cannot be too small 
since we expect the simulation to be finished in reasonable 
time. Nevertheless sometimes it is necessary to stabilize its 
output. One tested method is by keeping certain regions of 
the lattice with value of superconducting order parameter set 
to zero as it is depicted in Fig. 8. In order to confirm the 
stability of solution is the addition of noise to the system. 
After certain time the system should recover all values of 
A z (x,y) from before the noise addition. In this way we can 
become more sure that we have obtained the stable numerical 
solution. This is important since the space of initial probe 
as A z (x,y), i/r(x,y) functions is infinite. In general from the 
observation we have noticed that the probe function should 
have bigger monotonicity change than the expected numerical 
solution. We have tested the relaxation algorithm for the case 
of superconducting rectangular for d-wave superconductor 
in ab-plane. Topology of solutions of GL(v 2 — y 2 ) obtained 
by the relaxation method was in accordance with solutions 
obtained by different methods as described by [?]. 

Further perspectives 

The relaxation method used for study of superconducting 
mesoscopic structures is very stable even in the cases of more 
complex GL functionals as given in [?]. This method have the 


capacity to model very complex mesoscopic superconduct¬ 
ing structures. The parallelization of the relaxation method 
should be introduced. Because of its simplicity we have the 
reasons to believe that relaxation method could be the core for 
building universal platform capable of modeling many types 
of superconducting devices and various mesoscopic structures. 
Various results obtain in [?] needs to be further confirmed by 
the relaxation method. 
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Figure 9. Schematic description of relaxation algorithm solving GL/TDGL equations. 























































































